The Care and Feeding of 2D Data using Shapely¶
Overview and Examples¶
Shapely is a Python package that can be used to create and manipulate 2D Spatial data. This talk will just be an overview of some of the most useful functions available through the Shapely library. For more information, check out http://toblerity.org/shapely/manual.html¶
To get started, let's look at some potential applications of 2D data. For example, one could reimagine an implementation of Conway's game of life that is free of any grid.¶
In [192]:
from IPython.core.display import Image
Image(filename='files/conway.png')
Out[192]:
I use Shapely regularly when dealing with geospatial datasets, mostly in a geospatial context. Below is a map I created by layering data about different types of alcohol drinking habits. Can you guess what the darker patches in the southern and coastal parts of the state indicate?¶
In [15]:
from IPython.core.display import Image
Image(filename='files/GA.png')
Out[15]:
The example below is courtesy of Mike Higgins. He used our mapping implementation to draw a picture of Luna the Cat near the equator. In this case, he mapped latitiude and longitudes and then used a radius measure to draw circles around those points, forming a grid.¶
In [13]:
Image(filename='files/luna.jpg')
Out[13]:
Drawing Basic Shapes¶
In [24]:
import os
import numpy as np
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import pysal
import pylab as pl
import shapely
import shapely.geometry
from shapely.geometry import shape
from shapely.geometry import Point, LineString, Polygon, MultiPolygon
from shapely.wkt import dumps, loads
#from shapely.wkb import dumps, loads
The Code below will help us draw our shapes on the screen using matplotlib. (ADD SOURCING and ADDITIONAL COMMENTS)¶
In [292]:
import random
def get_random_color():
r = lambda: random.randint(0,255)
return('#%02X%02X%02X' % (r(),r(),r()))
def plot(shapelyGeometries, color_dict={"fill":"#AADDCC", "line":"#666666","hole_fill":"#ffffff", "hole_line":"#999999" }):
'Plot shapelyGeometries'
figure = pl.figure(num=None, figsize=(10, 10), dpi=180)
axes = pl.axes()
axes.set_aspect('equal', 'datalim')
axes.xaxis.set_visible(False)
axes.yaxis.set_visible(False)
draw(shapelyGeometries, color_dict)
def draw(gs, color_dict):
'Draw shapelyGeometries'
# Handle single and multiple geometries
try:
gs = iter(gs)
except TypeError:
gs = [gs]
# For each shapelyGeometry,
for g in gs:
gType = g.geom_type
if gType.startswith('Multi') or gType == 'GeometryCollection':
draw(g.geoms, color_dict)
else:
draw_(g, color_dict)
def draw_(g, color_dict):
'Draw a shapelyGeometry; thanks to Sean Gilles'
gType = g.geom_type
if gType == 'Point':
pl.plot(g.x, g.y, 'k,')
elif gType == 'LineString':
x, y = g.xy
pl.plot(x, y, 'b-')
elif gType == 'Polygon':
#can draw parts as multiple colors
if not color_dict:
color_dict={"fill":get_random_color(),
"line":"#666666",
"hole_fill":"#FFFFFF",
"hole_line":"#999999" }
x, y = g.exterior.xy
pl.fill(x, y, color=color_dict["fill"], aa=True)
pl.plot(x, y, color=color_dict["line"], aa=True, lw=1.0)
for hole in g.interiors:
x, y = hole.xy
pl.fill(x, y, color=color_dict["hole_fill"], aa=True)
pl.plot(x, y, color=color_dict["hole_line"], aa=True, lw=1.0)
Let's draw a point.¶
In [237]:
pt = Point(0, 0)
print pt
plot(pt)
Okay, that is not much to look at. ENHANCE!!¶
In [281]:
circle = pt.buffer(.5)
plot([circle, pt])
In [278]:
print circle
In [290]:
small_circle = pt.buffer(.25)
plot([circle, small_circle, pt], color_dict=None)
You can also use one polygon to make a hole in the other.¶
In [283]:
donut = circle.difference(small_circle)
plot(donut)
Let's draw a line. It's essentially made of two points.¶
In [242]:
lineString1 = LineString(((-3, -5), (3, 5)))
print lineString1
plot(lineString1)
How about a super crazy line? Oh yes, this is super crazy.
In [243]:
lineString2 = LineString(((-3, -5), (3, 5), (-3, 5), (3, -5)))
print lineString2
plot(lineString2)
And now we can add polygons to our linestrings to make a sad bunny.
In [244]:
pt1 = Point(-3, 6.5)
pt2 = Point(3, 6.5)
left = pt1.buffer(.5)
right = pt2.buffer(.5)
poly = shapely.geometry.Polygon(
((-3, 5), (3, 5), (0, 0), (-3, 5)))
plot([lineString2, left, right, poly], {"fill":"#000000", "line":"#AA8888","hole_fill":"#ffffff", "hole_line":"#999999" })
Let's look at a Polygon with holes.
In [245]:
polygon1 = shapely.geometry.Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 3), (3, 3), (3, 2), (2, 2)), # Hole
((3, 3), (3, 4), (4, 4), (4, 3), (3, 3)), # Hole
]
)
print polygon1
plot(polygon1)
You can use union() to create more complex shapes.
In [293]:
circle = Point(5,2.5).buffer(2)
plot([circle, polygon1], color_dict=None)
Now, let's look at some more complicated shapes.
In [247]:
new_polygon = polygon1.union(circle)
plot(new_polygon)
print new_polygon.type
In [294]:
circle2 = Point(7,7).buffer(2)
plot([circle2, polygon1], color_dict=None)
In [249]:
new_multipolygon = polygon1.union(circle2)
plot(new_multipolygon)
print new_multipolygon.type
print "CIRCLE AREA: %f" % (circle2.area)
print "POLYGON AREA: %f" % (polygon1.area)
print "MULTIPOLYGON AREA: %f" % (new_multipolygon.area)
Getting Shapes from a File¶
One other way to get shape data is to get it from a file. Below is an example of a function that creates a dictionary of shapes keyed by an id column
In [250]:
import fiona
def get_from_fiona(filename, id_attr):
from shapely.geometry import shape
data = {}
c = fiona.open(filename, 'r')
for shape_data in c:
name = shape_data['properties'][id_attr]
parsed_shape = shape(shape_data['geometry'])
if (name in data):
print ("DUPLICATE: " + str(name))
dup_shape = data[name]
new_shape = merge_duplicates(shapely.wkt.dumps(dup_shape), shapely.wkt.dumps(parsed_shape))
parsed_shape = new_shape
data[name] = parsed_shape
#print name + " : " +str(parsed_shape)
return data
Just a quick note, you can get your own shape data from the US Census here: https://www.census.gov/cgi-bin/geo/shapefiles2013/main
We're going to look at state shapes from the US Census. Just to give you an idea of what is in the file, here is a glance at the 'properties' of a few records.
In [251]:
from shapely.wkt import dumps, loads
state_dict = get_from_fiona("files/tl_2013_us_state/tl_2013_us_state.shp", "NAME")
#print(shapely.wkt.dumps(state_dict["Ohio"]))
ohio = state_dict["New Jersey"]
print ohio.convex_hull
In [252]:
plot(state_dict["Pennsylvania"])
print
In [297]:
plot([state_dict["Pennsylvania"],
state_dict["Ohio"],
state_dict["West Virginia"],
state_dict["Maryland"],
state_dict["Kentucky"],
state_dict["Indiana"],
state_dict["Michigan"]], color_dict=None)
print
You might notice that Michigan and Maryland look a little/lot wrong. This is because the USGS sometimes includes water features that are counted as part of a state's territory as part of their borders. Maryland has the Chesapeake as part of its shape and Michigan includes a great deal of the Great Lakes. Below is a more familiar picture of what Michigan looks like.
In [254]:
Image(filename='files/MI.png')
Out[254]:
Simplifying Shapes¶
It's also handy to simplify shapes. This can be a net good for speed, especially where details aren't as pertinent. However, there are drawbacks in the form of loss of precision and potential introduction of unsound polygons. Shapely's simplification methodology is based on the Ramer-Douglas-Peucker algorithm. For more information see here: http://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
"Douglas-Peucker animated" by Mysid - Own work; self-made in Inkscape and Gimp. Based on File:Douglas Peucker.png by de:User:Leupold.. Licensed under CC BY-SA 3.0 via Wikimedia Commons.
In [255]:
texas_simple = texas.simplify(.001)
plot(texas_simple)
print "TOTAL POINTS: %d" % (points(texas_simple))
In [256]:
texas_simpler = texas.simplify(.05)
plot(texas_simpler)
print "TOTAL POINTS: %d" % (points(texas_simpler))
In [257]:
texas_simplest = texas.simplify(.1)
plot(texas_simplest)
print "TOTAL POINTS: %d" % (points(texas_simplest))
In [258]:
not_texas = texas.simplify(1.1)
plot(not_texas)
print "TOTAL POINTS: %d" % (points(not_texas))
Let's look at some further operations that can help us simplify and understand the shapes we are generating.
In [315]:
plot(texas.convex_hull)
In [312]:
plot(texas.envelope)
In [314]:
print(texas.bounds)
In [322]:
print(texas.centroid)
center_circle1 = texas.centroid.buffer(.2)
center_circle2 = texas.centroid.buffer(.5)
center_circle3 = texas.centroid.buffer(1)
plot([texas, center_circle3, center_circle2, center_circle1, texas.centroid], color_dict=None)
More Polygon Operations¶
In [303]:
from shapely.ops import cascaded_union, unary_union
buffered_points_list = []
x_ys = range(5);
for x in x_ys:
buffered_points_list.append(Point(x,x).buffer(1))
union_shape = cascaded_union(buffered_points_list)
plot(buffered_points_list, color_dict=None)
plot(union_shape)
Intersections are great for showing the degree of overlap as well as creating new shapes.
In [300]:
intersection_shape = buffered_points_list[0].intersection(buffered_points_list[1])
plot([buffered_points_list[0], buffered_points_list[1], intersection_shape], color_dict=None)
Difference is the inverse of intersection and shows the parts of a shape that do not overlap. Unlike the intersection function, this operation is not symmetrical
In [301]:
difference_shape1 = buffered_points_list[0].difference(buffered_points_list[1])
difference_shape2 = buffered_points_list[1].difference(buffered_points_list[0])
plot([difference_shape1,difference_shape2], color_dict=None)
Here are a few examples of shape tests just to show how subtle differences in test definitions can affect results.
In [329]:
print "Do the full circle shapes intersect?"
print(buffered_points_list[0].intersects(buffered_points_list[1]))
print "Do the difference shapes intersect?"
print(difference_shape1.intersects(difference_shape2))
print "Are the difference shapes overlapping?"
print(difference_shape1.overlaps(difference_shape2))
Finding Broken Shapes and Fixing Them¶
First, there are self-intersections
In [261]:
sad_polygon1 = Polygon(((-3, -5),(-3,0),(-5,0),(-3,0),(-3, 5),(3, 5), (3, -5), (-3, -5)))
plot(sad_polygon1)
print explain_validity(sad_polygon1)
In [262]:
print "SAD POLYGON AREA: %f" % (sad_polygon1.area)
fixed_polygon1 = sad_polygon1.buffer(.03)
plot(fixed_polygon1)
print explain_validity(fixed_polygon1)
print "HAPPY POLYGON AREA: %f" % (fixed_polygon1.area)
Let's look at another type of self-intersection, the bowtie. It can be more tricky than dealing with the extra line.
In [263]:
from shapely.validation import explain_validity
sad_polygon2 = Polygon(((-3, -5), (3, 5), (-3, 5), (3, -5), (-3, -5)))
plot(sad_polygon2)
print explain_validity(sad_polygon2)
Using the buffer method for cleanup can lead to a deletion of part of the unsound shape. Notice that the area of this shape is calculated as 0.
In [264]:
print sad_polygon2.area
fixed_polygon2 = sad_polygon2.buffer(0)
plot(fixed_polygon2)
print explain_validity(fixed_polygon2)
print fixed_polygon2.area
Other solutions also fail due to the fact that the shape is too unsound for unions.
In [302]:
from shapely.validation import explain_validity
print explain_validity(sad_polygon2)
pt = Point(0,0)
band_aid = pt.buffer(.2)
plot([sad_polygon2, band_aid], color_dict=None)
try:
trythis = sad_polygon2.union(band_aid)
except(shapely.geos.TopologicalError):
print "ERROR"
Here is a quick look at the stack trace without the exception handling for your edification.
In [266]:
---------------------------------------------------------------------------
TopologicalError Traceback (most recent call last)
<ipython-input-157-cfbdb28f0e7f> in <module>()
4 band_aid = pt.buffer(.2)
5 plot([sad_polygon2, band_aid])
----> 6 trythis = sad_polygon2.union(band_aid)
7
8 try:
/Users/alisonalvarez/anaconda/lib/python2.7/site-packages/shapely/geometry/base.pyc in union(self, other)
489 def union(self, other):
490 """Returns the union of the geometries (Shapely geometry)"""
--> 491 return geom_factory(self.impl['union'](self, other))
492
493 # Unary predicates
/Users/alisonalvarez/anaconda/lib/python2.7/site-packages/shapely/topology.pyc in __call__(self, this, other, *args)
45 if not this.is_valid:
46 raise TopologicalError(
---> 47 "The operation '%s' produced a null geometry. Likely cause is invalidity of the geometry %s" % (self.fn.__name__, repr(this)))
48 elif not other.is_valid:
49 raise TopologicalError(
TopologicalError: The operation 'GEOSUnion_r' produced a null geometry. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x110846350>
The best way to deal with this kind of unsound shape is to not make it in the first place. In the case where this isn't really an option (such as with simplification) it's possible to clean up and simplify component shapes to prevent this. Complexity should be tamed.
In []:
Let's take a look at some unsound shapes that are a result of holes in the wrong place.
In [267]:
sad_polygon3 = shapely.geometry.Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 3), (3, 3), (3, 2), (2, 2)), # Hole
((6, 6), (6, 7), (7, 7), (7, 6), (6, 6)), # Hole
]
)
plot(sad_polygon3)
print explain_validity(sad_polygon3)
In [268]:
print "SAD POLYGON AREA: %f" % (sad_polygon3.area)
fixed_polygon3 = sad_polygon3.buffer(0)
plot(fixed_polygon3)
print explain_validity(fixed_polygon3)
print "FIXED POLYGON AREA: %f" % (fixed_polygon3.area)
print fixed_polygon3
In [269]:
sad_polygon4 = shapely.geometry.Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 4), (4, 4), (4, 2), (2, 2)), # Hole
((3, 3), (3, 3.5), (3.5, 3.5), (3.5, 3), (3, 3)), # Hole
]
)
plot(sad_polygon4)
print explain_validity(sad_polygon4)
In [270]:
print "SAD POLYGON AREA: %f" % (sad_polygon4.area)
fixed_polygon4 = sad_polygon4.buffer(0)
plot(fixed_polygon4)
print explain_validity(fixed_polygon4)
print "FIXED POLYGON AREA: %f" % (fixed_polygon4.area)
print fixed_polygon4
In [271]:
sad_polygon5 = shapely.geometry.Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 3.5), (3, 3.5), (3, 2), (2, 2)), # Hole
((3, 3), (3, 4), (4, 4), (4, 3), (3, 3)), # Hole
]
)
print explain_validity(sad_polygon5)
print sad_polygon5.is_valid
plot(sad_polygon5)
In [272]:
print "SAD POLYGON AREA: %f" % (sad_polygon5.area)
fixed_polygon5 = sad_polygon5.buffer(0)
plot(fixed_polygon5)
print explain_validity(fixed_polygon5)
print "FIXED POLYGON AREA: %f" % (fixed_polygon5.area)
print fixed_polygon5
In [273]:
sad_polygon6 = shapely.geometry.Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 2), (3, 3), (3, 3), (2, 2)), # Hole
]
)
plot(sad_polygon6)
print explain_validity(sad_polygon6)
In [274]:
def remove_bad_shapes(shape):
if shape.type == 'Multipolygon':
shape_list = shape.geoms
elif shape.type == 'Polygon':
shape_list = [shape]
else:
return None
polygon_list = []
for g in shape_list:
print g
print explain_validity(g)
hole_list = []
for hole in g.interiors:
hole_shape = Polygon(hole)
if "Too few points" not in explain_validity(hole_shape):
hole_list.append(hole)
else:
print "Removed hole: "
print "\t" + str(hole)
if "Too few points" not in explain_validity(Polygon(g.exterior.coords)):
print explain_validity(Polygon(g.exterior.coords))
polygon_list.append(Polygon(g.exterior.coords, hole_list))
else:
print "Removed: "
print "\t" + str(g.exterior.coords)
return MultiPolygon(polygon_list)
In [275]:
print "SAD POLYGON AREA: %f" % (sad_polygon6.area)
fixed_polygon6 = remove_bad_shapes(sad_polygon6)
plot(fixed_polygon6)
print explain_validity(fixed_polygon6)
print "FIXED POLYGON AREA: %f" % (fixed_polygon6.area)
A Closer Look at Cleaning Up Shapes¶
At Rhiza we do a lot of work generating novel shapes from many components of varying quality.
In [204]:
Image(filename='files/unsimplified_clayton.png')
Out[204]:
This shape is the result of removing exterior and interior rings that are below a certain size threshold.
In [203]:
Image(filename='files/clayton_simple.png')
Out[203]:
Here is another example from the union of zip codes.
In [223]:
Image(filename='files/separate_zips.png')
Out[223]:
In [222]:
Image(filename='files/not_cleaned.png')
Out[222]:
In [224]:
Image(filename='files/cleaned_up.png')
Out[224]:
In []: